/* Function for rapid calculation of 1D density of points at those points.
Assumes a Gaussian kernel.
Requires an array of 1D locations (X[i]) and a bandwidth (bw)
*/

#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
#include <Rdefines.h>
#include <math.h>

double Pi = 3.141593;
//--------------------------------------//
// Function declarations //

double norm (double x);
SEXP dens (SEXP R_X, SEXP R_n);


//--------------------------------------//
//  Function definitions  //

double norm (double x){
	double bw=1;
	return exp(-pow(x,2)/(2*pow(bw,2)))/(bw*sqrt(2*Pi));
}

SEXP dens (SEXP R_X, SEXP R_n){
	int n = INTEGER(coerceVector(R_n, INTSXP))[0]; //grab vector length
	PROTECT(R_X=coerceVector(R_X, REALSXP)); //digest R_X
	double *X, *D, d; //variables
	SEXP out; PROTECT(out = allocVector(REALSXP,n)); //a place to put the answer
	D = REAL(out); //pointer to output
	X = REAL(R_X); //pointer to vector

	for (int ii=0; ii<n; ii++){
		d = 0;
		for (int jj=0; jj<n; jj++){
			d+= norm(X[ii]-X[jj]);
		}
		D[ii] = d;
	}
	UNPROTECT(2);
	return(out);
	//Free(X);
	//Free(D);
	
	//Free(out);
}


